Processing sonic waveform measurements

ABSTRACT

Methods for determining the best value for at least one slowness-related parameter that has been determined in a number of ways is disclosed. Sonic logging data input in the methods are processed to determine multiple values of at least one slowness-related parameter using slowness/time coherence (STC) processing methods. The error of each determined parameter is determined and the determined errors used in selecting a representative parameter value from the multiple determined slowness-related parameter values.

FIELD OF THE INVENTION

The present invention relates to methods for processing sonic waveformmeasurements, particularly sonic waveform measurements made for thepurpose of characterising properties of underground formations. Theinvention in particular relates to methods for determining the bestvalue for a parameter that has been determined in a number of differentways.

RELATED ART

It has been know for some time that it is possible to determineproperties of underground formations using measurements ofacoustic/sonic waves that have passed through the formations. The basictechnique comprises placing a tool comprising a spaced sonic source andreceiver in the borehole and using the source to generate sonic waveswhich pass through the formation around the borehole and are detected atthe receiver. Sonic waves can travel through rock formations inessentially two forms: body waves and surface waves. There are two typesof body waves that travel in rock: compressional and shear.Compressional waves, or P-waves, are waves of compression and expansionand are created when rock formation is sharply compressed. Withcompressional waves, small particle vibrations occur in the ,samedirection the wave is travelling. Shear waves, or S-waves, are waves ofshearing action as would occur when a body is struck from the side. Inthis case, rock particle motion is perpendicular to the direction ofwave propagation. The surface waves are found in a borehole environmentas complicated borehole guided waves which come from reflections of thesource waves reverberating in the borehole. The most common form orborehole-guided, surface wave is the Stoneley wave. FIG. 1 shows aseries of sonic waveforms such as would be recorded in a borehole from amonopole (omnidirectional) source with the first arrivals of thecompressional (P), shear (S) and Stoneley (St) waves at the receivermarked. In situations where dipole (directional) sources and receiversare used, an additional shear/flexural wave propagates along theborehole and is caused by the flexing action of the borehole in responseto the dipole signal from the source The flexural wave typically travelsat about the same speed as the shear wave, slower than the compressionalwave. (It is to be noted that sonic waves will also travel through thefluid in the borehole and along the tool itself. With no interactionwith the formation, these waves carry no useful information and run onwireline or coiled tubing or the like, or alternatively can be a loggingwhile drilling tool located in a drill string being used to drill theborehole.

The speeds at which these waves travel through the rock are controlledby rock mechanical properties such as density and elastic dynamicconstants, and other formation properties such as amount and type offluid present in the rock, the makeup of the rock grains and the degreeof intergrain cementation. Thus by measuring the speed of sonic wavepropagation in a borehole, it is possible to characterise thesurrounding formations by parameters relating these properties. Thespeed or velocity of a sonic wave is often expressed in terms of1/velocity and is called “slowness”. Since the tools used to make sonicmeasurements in boreholes are of fixed length, the difference in time(ΔT) taken for a sonic wave to travel between two points on the tool isdirectly related to the speed/slowness of the wave in the formation.

An example of a tool for use in a borehole for sonic measurements is theDSI tool of Schlumberger which is shown schematically in FIG. 2. The DSItool comprises a transmitter section 10 having a pair of (upper andlower) dipole sources 12 arranged orthogonally in the radial plane and amonopole source 14. A sonic isolation joint 16 connects the transmittersection 10 to a receiver section 18 which contains an array of eightspaced receiver stations, each containing two hydrophone pairs, oneoriented in line with one of the dipole sources, the other with theorthogonal source. An electronics cartridge 20 is connected at the topof the receiver section 18 and allows communication between the tool anda control unit 22 located at the surface via an electric cable 24. Withsuch a tool it is possible to make both monopole and dipolemeasurements. The DSI tool has several data acquisition operating modes,any of which may be combined to acquire (digitised) waveforms. The modesare: upper and lower dipole modes (UDP, LDP)—waveforms recorded fromreceiver pairs aligned with the respective dipole source used togenerate the signal; crossed dipole mode—waveforms recorded from eachreceiver pair for firings of the in-line and crossed dipole source;Stoneley mode—monopole waveforms from low frequency firing of themonopole source; P and S mode (P&S)—monopole waveforms from highfrequency firing of the monpole transmitter; and first motionmode—monopole threshold crossing data from high frequency firing of themonopole source.

One way to determine compressional, shear and Stoneley slownesses fromthese measurements is to use slowness-time-coherence (STC) processing.STC processing is a full waveform analysis technique which aims to findall propagating waves in the composite waveform. The processing adopts asemblance algorithm to detect arrivals that are coherent across thearray of receivers and estimates their slowness. The basic algorithmadvances a fixed-length time window across the waveforms in small,overlapping steps through a range of potential arrival times. For eachtime position, the window position is moved out linearly in time, acrossthe array of receiver waveforms, beginning with a moveout correspondingto the fastest wave expected and stepping to the slowest wave expected.For each moveout, a coherence function is computed to measure thesimilarity of the waves within the window. When the window time and themoveout correspond to the arrival time and slowness of a particularcomponent, the waveforms within the window are almost identical,yielding a high value of coherence. In this way, the set of waveformsfrom the array is examined over a range of possible arrival times andslownesses for wave components. STC processing produces coherence(semblance) contour plots in the slowness/arrival time plane. Regions oflarge coherence correspond to particular arrivals in the waveforms. Theslowness and arrival time at each coherence peak are compared with thepropagation characteristics expected of the arrivals being sought andthe ones that best agree with these characteristics are retained.Classifying the arrivals in this manner produces a continuous log ofslowness versus depth. For dispersive waves, the STC processing ismodified to take into account the effect of frequency. As the output ofSTC processing is a coherence plot, the coherence of each arrival can beused as a quality indicator, higher values implying greater measurementrepeatability. When processing dipole waveforms, one of the coherencepeak will correspond to the flexural mode but with a slowness that isalways greater (slower) than the true shear slowness. A precomputedcorrection is used to remove this bias.

To compensate for variations in measurements due to the borehole ratherthan due to the formation a series of measurements are made across aninterval in which the formation properties are expected to vary little,if at all. In its simplest form, the interval corresponds to the extentof the receiver array, and the waveforms at each receiver stationmeasured for a given firing of a source (“receiver array” or “receivermode” or “Rec.”). In simple STC processing, all receiver stations areconsidered. In multishot STC processing (MSTC), sub-arrays of receiverstations within the receiver array are considered, for example asub-array of five receiver stations in a receiver array of eightreceiver stations (other numbers or receiver stations in the sub-arraycan be used depending on requirements). In this case, the same formationinterval corresponding to the extent of a five receiver stationsub-array can be measured several times as the tool is logged throughthe borehole, the five stations making up the sub-array being selectedat each source firing to measure the same formation interval. Anotherapproach, known as “transmitter mode” or “pseudo-transmitter array”(“Tra.”) takes waveforms from sequential source firings as thetransmitter passes along the interval to be measured. In order tocompensate for the movement of the tool between measurements, aneffectively stationary receiver station or sub-array must be used. Thiscan be achieved by changing the receiver station considered so that itsposition in the borehole is effectively stationary as the transmitter ismoved through the interval. Borehole compensation (“BHC”) can beachieved for P and S mode results by processing receiver array andpseudo-transmitter array waveforms and averaging the results.

Thus it will be appreciated that, with the different acquisition modesof a tool such as the DSI, and the different processing modes available,it is possible to obtain multiple determinations of a slowness or ΔT ina given interval of borehole. For example, it is possible to acquirewaveforms for shear slowness determination using two dipole modes (upperdipole and lower dipole), and one monopole mode (P and S mode), and toprocess each measurement in receiver mode, transmitter mode and boreholecompensated form resulting in a potential nine separate determinationsof shear slowness for that interval, each of which can give a differentresult. The problem is therefore to determine which slowness estimatecan be considered to give the best indication of the shear slowness ofthe formation in that interval.

DISCLOSURE OF INVENTION

The present invention provides a method of determining the sonicslowness of an underground formation from sonic measurements,comprising: (i) obtaining sonic waveforms in the underground formation;(ii) determining, from the sonic waveforms, multiple values of at leastone parameter related to the sonic slowness of the formation togetherwith an estimate of the error in each value; and (iii) using theestimate of the error in each value to select a parameter value relatedto slowness as representative of the formation.

The sonic waveforms are preferably obtained by logging an interval of aborehole which runs through the formation with a tool which outputssonic waveform measurements. The tool can be run on wireline or coiledtubing or the like, or alternatively can be a logging while drillingtool located in a drill string being used to drill the borehole.

The multiple values of the parameter related to sonic slowness caninclude multiple determinations of monopole and/or dipole compressionaland shear (including flexural) slowness, and Stoneley slowness for theformation. Where the tool used to obtain the waveforms comprise an arraytool, the multiple determinations can include transmitter and receivermode measurements and borehole compensated measurments. The processingtechnique used is preferably a slowness-time-coherence technique.

The processing technique can include amongst its inputs, zoninginformation derived from the waveform measurements and indicatinggeneral features of the formation type being measured. The zoninginformation can be obtained from a basic compressional slownessestimation, typically based on digital first arrival determination, andcan include broad formation slowness classifications such as fast, slow,very slow and extremely slow. Such broad classifications can be based onpredetermined cross plots of the ratio of compressional and shearslownesses against measured slowness for know lithologies. Other zoninginformation can be the presence of closely spaced bed boundaries (thinbeds).

Other inputs to the processing technique can include parameters relatingto the borehole or well, such as hole diameter, mud type andpredetermined formation features.

The specific processing technique applied can vary according to thezoning or parameter inputs.

For example, where zoning information shows relatively thick beds, fullarray STC-type processing (including dispersive processing) can beapplied; if zoning information shows thin beds, high resolution,sub-array multishot STC-type processing can be applied. The processingalso preferably includes tracking of slowness measurements along theinterval so as to allow individual slowness measurements to beassociated with changes in particular components of the waveform(P-waves, S-waves, etc.). The tracking can also make use of the zoninginformation indicating where major changes occur, and delineatinghomogenous beds, each with associated semblance error bars.

The output of the processing step is a series of slowness estimates andthe tracking can also include the estimation of error in the value ofslowness. An estimate of error can be provided by the statisticalsemblance processing of the slowness measurements. The total error is acombination of the semblance error and the uncertainty in trackingdetermination. The final step selects the slowness with the minimumerror (possibly modified by other predetermined selection rules) andoutputs this as the slowness for the interval. For example, the meanerror for the interval can be used as the basis for selection. Also, alevel by level determination of the variation of the actual error fromthe mean within the interval can be provided as a further output

BRIEF DESCRIPTION OF DRAWINGS

FIGS. 1 shows a plot of sonic waveforms in time indicating arrivals ofvarious components;

FIG. 2 shows a prior art sonic logging tool;

FIG. 3 shows a flow diagram outlining the main features of a methodaccording to one embodiment of the invention;

FIG. 4 shows a simple compressional slowness log;

FIG. 5 shows a Vp/Vs vs compressional slowness cross plot; and

FIG. 6 shows a slowness log with multiple slowness determinations.

DETAILED DESCRIPTION

Referring now to FIG. 3, the method summarised therein takes as itsinputs the raw waveforms 100 obtained from a borehole logging tool suchas the DSI tool from Schlumberger, zoning information 110 and processingparameters 120. The raw waveforms 100 are essentially digital signalscomprising the receiver station output in time for a given acquisitionmode of the tool (cf. FIG. 1). To obtain the zoning information 110, theraw waveforms 100 are pre-processed using a digital first arrivaldetection (DFAD) method such as is described in more detail inInternational Patent Application No. WO97/28464 to derive a crudecompressional slowness (DTCO) log vs depth 102 an example of which isshown in FIG. 4. The DTCO log information is analysed 104 to determinethe general zones present in the interval of interest. This is achievedby applying predetermined thresholds to the DTCO log and squaring theoutput so as to broadly classify the log into zones of fast, slow, veryslow and extremely slow formation with sharp transitions between thezones (line DTCOsq in FIG. 4), and is called “macro-zoning”. Where thezones are changing frequently over relatively short distances (thinbeds, for example <2 ft thickness), the log is also indicated as having“micro-zoning”. The macro- and micro-zoning information is output to aΔT computation process 130 (micro-zoning information can be used toindicate bed boundaries and to select multishot STC processing which isdescribed in more detail below) and (for macro-zoning information)to aparameter selection process 112.

The parameter selection process 112 determines from physical wellparameters 114 a series of processing parameters (see examples given inthe table below) and provides these as a processing parameter input 120to the ΔT computation process 130.

Well Parameter Tool Type Formation Type Select from automatic usingzoning input or manual selection based on VpVs-DTCO crossplot overlaidwith formation type regions. Hole Diameter Bit size, casing ID orcaliper Borehole Status Open Hole or Cased Hole Processing ParameterDispersion Curve Slowness dispersion curve for the mode to be evaluated.Integration Time Window Window size for coherence calculation, fewcycles of lowest frequency component Slowness Lower Limit Start, stopand increment along Slowness Upper Limit slowness axis of S/T plane forSTC Slowness Step Parameter processing. Time Lower Limit Start, stop andincrement along Time Upper Limit time axis of S/T plane for STC TimeStep Parameter processing. Search Band Width Define limits of S/T planeto be Search Band Offset processed. Slowness Width Define mask for peaksearch. Time Width Semblance Threshold Minimum semblance value forpeaks. Filter Length, Upper Limit, Lower Processing frequency band.Limit Mud Global Parameters Mud type: Brine, WBM, OBM, Slow OBM; orcalculate from mud type (brine, emulsion, invert emulsion) salinity,density, temperature.

It will be appreciated that not all of these parameters might be needed,or that default values might be acceptable in some cases. The particularparameters, how they are derived and their significance will depend onthe, particular processing scheme used. There are other parameters whichare optional such as optimal frequency selection for processing such asis described in more detail in U.S. Pat. No. 5,587,966.

The parameter “Formation Type” can have two main methods of selection,both of which rely on the zoning step 104. In the first, the DTCO logdata for a given level is plotted on a scale which is banded intopredefined formation types. In FIG. 5, DTCO data from two depths areindicated by + and x. For the purpose of illustration, these are crossplotted with Vp/Vs and the predefined formation types: fast (F),intermediate (I), slow (S), very slow (VS) and extremely slow (XS)superimposed on the plot The predefined formation types are derived fromother log data obtained in a number of locations. The data from depth +fall in a fast formation type while those from depth x fall in a slowformation type. Thus different parameters will be selected forprocessing the data from these two depths. It will be appreciated that,normally, there will be no Vp/Vs data available and it is only positionon the DTCO axis which is used to decide formation type. An alternativemethod is to output DTCOsq directly as the formation type.

The ΔT computation process 130 comprises two main elements, anintegrated slowness determination process 132 and a tracking process134. The integrated slowness determination process 132 provides aSlowness/Time Coherence (STC) methods for analysing the waveforms. Thesemethods can include one or more of STC, multishot STC (MSTC), fast STC(FSTC) and dispersive STC (DSTC). These different methods are summarisedbelow. In the preferred example, the major processing method utilised isDSTC.

Slowness Time Coherence (STC) Processing

STC processing employs full waveform analysis to find all propagatingwaves in an array of sonic waveforms. The processing computes a peakvector that identifies all of the arrivals in the waveform data from asonic array tool such as the DSI. The peak vector consists of sevenelements associated with a peak coherence value in the slowness/timeplane. Two main steps are carried out in the STC technique:

1. The coherence σ(T,S) is calculated for all reasonable values of Timeand Slowness in the S/T plane. σ(T,S) serves as a measure for whetherthe data includes an arrival at time T and slowness S.

2. Then, the surface σ(T,S) is searched for local maxima with a peakfinding process.

The coherence measure is calculated by:${\sigma \left( {T,S} \right)} = \sqrt{\frac{E_{c}\left( {T,S} \right)}{{MxE}_{i}\left( {T,S} \right)}}$

Where E_(c)(T,S) is the coherent energy of the normalised waveforms(square of the summed waveform values), E_(i)(T,S) is the total energyof the normalised waveforms (sum of the squared waveform values), and Mis the number of waveforms, i.e. receivers (receiver stations).

E_(c)(T S) and E_(i)(T,S) are calculated over a specified time window,usually selected to span several cycles of the main frequency componentexpected in the waveforms. In E_(c)(T,S), the waveforms are stacked intime assuming the arrivals move out linearly across the receiver array.For any number of waveforms, the coherence measure σ(T,S) will bebetween 0 and 1. Values approaching 1 indicate the presence of a highlycoherent arrival at time T and slowness S.

The coherence σ(T,S) is calculated within a restricted search band ofthe S/T plane defined by the slowness lower and upper limits and timelower and upper limits. For a given slowness S between the lower andupper slowness limits, the time (T) in the search band must obey thefollowing rule:

S·Z0+T _(off) −T _(band)/2≦T≦S·Z0+T _(off) +T _(band)/2

where Z0 is the T/R spacing of the tool, T_(off) is the time offset ofthe search band, and T_(band) is the width along time of the searchband.

This rule describes a band of times moving diagonally across the S/Tplane from the lower, left-hand corner to the upper, right-hand corner.For T_(off) equal to 0, this band is centred on the slowness-time line(T=S·Z0).

The point (T_(pk), S_(pk)) is defined as the peak coherence in the S/Tplane search ask if it satisfies the following:

1. σ(T_(pk), S_(pk))>0.25 (Threshold value) A rule of thumb is thatvalues greater than 2/M indicate the presence of an arrival. If 8waveforms are being input, the value of 0.25 can be assumed.

2. σ(T_(pk), S_(pk)) is a local maximum within the given peak (specifiedby the slowness width an time width):

σ(T _(pk) , S _(pk))>p(T _(m) , S _(m)) for all T_(m) and S_(m)

 such that

T−T_(mask)/2≦T_(m)≦T+T_(mask)/2

S−S_(mask)/2≦S_(m)≦S+S_(mask)/2

where T_(mask) is the time width of the peak mask, and S_(mask) is theslowness width of the peak mask.

The following peak matrix of seven elements are output to seven arraysat every depth processed:

At peak coherence:

1. σ(T_(pk), S_(pk)) coherence value

2. S_(pk) slowness

3. T_(pk) time

4. E_(c)(T_(pk), S_(pk)) coherent energy in decibels

At peak energy:

5. σ(T_(e), S_(pk)) coherence

6. T_(e) time

7. E_(c)(T_(e), S_(pk)) coherent energy in decibels

The slowness is the same as for peak coherence.

Signal-to-noise ratio is estimated using the following algorithm:

1. The peak-to-peak amplitude in the signal window is determined at eachreceiver and the amplitudes are averaged to get an estimate of thesignal.

2. The peak-to-peak amplitude in the quite window is determined at eachreceiver and the amplitudes are averaged to get an estimate of thenoise.

3. The signal-to-noise ratio estimate is then simply calculated usingthe following equation:

SIN=20·Log 10 (signal/noise).

There are two main depth reference points in the STC processing (or DSTCprocessing, see below), waveform depth and computation depth. Waveformdepth is the depth assigned to all waveforms in the array at aparticular depth level and is assumed to be the mid point between thefired transmitter and the first receiver in the receiver array. Allwaveforms in the array are assigned this depth. The computation depth isthe depth assigned to the S/T plane and the peak matrix computed by oneSTC computation. When a receiver array mode is selected, this depth isthe mid-point of the receiver array used in the computation. Thus if alleight receivers are used, the depth is the mid-point between receivers 1and 8. If a sub-array of receivers 5-8 is used, the depth is themid-point between receivers 5 and 8. When a transmitter array mode isused, the depth is the mid-point of the transmitter positions formingthe pseudo-array.

Dispersive STC (DSTC) Processing

Certain modes of sonic waves are dispersive, for example the dipoleflexural shear mode. This means that different frequencies propagate atdifferent phase slownesses, waveshapes change as they propagate acrossthe array, and measurement of shear slowness requires a dispersionmodel. Processing the flexural mode with standard, non-dispersiveprocessing such as STC (described above) does not give formation shearslowness without the use of model-based corrections. Dispersive STC(DSTC) processes the flexural mode dispersively so no correction fordispersion is required. In STC, waveforms in the window arebackpropagated to the reference receiver along lines of constantslowness. This constant slowness is the same slowness as the moved-outwindow and the calculation can be done entirely in the time domain. Thesemblance of the backpropagated waveforms is calculated, high semblancevalues indicating that the backpropagated waveforms are nearly alikewhich, in turn, indicates that the backpropagating slowness is correct.DSTC performs backpropagation in a different way. The waveforms arebackpropagated by a flexural mode dispersion curve equal to that of thedata. The dispersion curve is pre-computed from the flexural mode modeland for speed the backpropagation is performed in the frequency domain.The shear slowness of the best-fitting curve is returned by DSTC. DSCTalso outputs an error bar for the uncertainty of the slownessdetermination due to noise. Further details of DSTC can be found in U.S.Pat. No. 5,278,805.

DSTC processing has the ability to be applied even to non-dispersivewaves such as monopole compressional or shear-. Since the first steprequired for DSTC processing is the calculation or selection of anappropriate dispersion curve, all that is required is a dispersion curvethat represents a non-dispersive wave, i.e. a flat “curve”.Consequently, DSTC is the preferred processing scheme for the presentinvention and will be indicated as such in the following description.However, it will be appreciated that in many cases, STC processing couldequally be applied.

Multishot STC (MSTC) Processing

MSTC attempts to improve the resolution of the DSTC full arraycalculation by processing results from sub-arrays and oversampling theinterval of interest. For example, in an eight receiver tool, it ispossible to improve resolution by forming smaller sub-arrays of, forexample, five receivers. While this reduces the amount of data availablefrom each sub-array, the fact that the tool has multiple opportunitiesto sample the region with different sub-arrays of the same size meansthat this can be compensated in the redundancy of the data. MSTC doesthis using the following procedure:

1. All of the sub-arrays which span the interval of interest areidentified.

2. Each sub-array is processed using DSTC (an S/T plane is generated foreach sub-array).

3. The S/T planes for the sub-arrays are shifted in time to account forthe difference in transmitter-receiver spacing.

4. The semblances from the various sub-arrays are averaged at everyvalue of time and slowness to provide and average S/T plane.

5. The average S/T plane is searched for slowness peaks which arereturned in a peak matrix.

The peak matrix output by MSTC contains three elements:

the slowness of the coherence peak,

the coherence of the coherence peak, and

the time of the coherence peak.

The slowness, time and coherence are simply the slowness, time andcoherence of each peak found in the average S/T plane.

Selection of MSTC as the processing is made when the zoning information110 indicates thin beds (typically <2 ft thickness). A single-shot DSTCcomputation by MSTC is the same as an STC computation made by DSTC.Further details of MSTC can be found in U.S. Pat. No. 4,809,236. Thegeneral processing methodology preferred for the present invention isdescribed in further detail below.

Fast STC (FSTC) Processing

Fast STC (FSTC) processing was developed for use with LWD sonicmeasurement techniques but is potentially an alternative to the standardSTC processing described above. In FSCT the raw waveforms are filteredto give a complex signal type with both real and imaginary parts.Because of the increase in data arising from this, the data are coarselyresampled to reduce the amount of data and only the real part of thecomplex signal data are saved for display purposes. The basic processingother than the filtering and sampling are as for STC. Further details ofFSTC can be found in U.S. Pat. No. 5,594,706.

Whichever processing method is used, the output is typically more thanone slowness for each depth, e.g. compressional slowness, shear slownessand Stoneley slowness. FIG. 6 shows a log with more than one slowness(e.g. compresional, shear, Stoneley) indicated for each level. It isgenerally desired to monitor the development of these properties alongthe well and so it is necessary to decide which slowness is which fromdepth to depth in the well. This can be difficult in cases where theslownesses are similar in value or cross each other, or are missingbecause of bad data at a given depth. Commonly used techniques forconnecting peaks in a continuous curve are called labelling and some usea logical decision tree, level by level. More sophisticated techniquesuse a weighted cost function over a short interval to select an optimumpath. The particularly preferred approach is to use a two-step process134, first joining the peaks that correspond to the same waveformarrival (“track search”), and second identifying the track with a name(“classification”).

To perform the track search, the sequence of measurements that aresufficiently continuous are enumerated and their probabilities computed.In order to know how to associate the different peaks with the differenttrack built in this way, all possible hypotheses are computed and thecorresponding probabilities derived. With this information, and byapplying Bayes theorem, the a posteori probability can be computed andthe hypothesis with the most probable value selected to build a track.FIG. 6 shows point on the log connected to form tracks. Classificationis achieved by forming different hypotheses based on a propagation modelto give all possible classifications of tracks built, and associating apriori probabilities with these hypotheses. Knowing the data and the apriori probability model, the a posteori probabilities of the differenthypotheses can be computed and the most probable hypothesis consideredto classify the tracks 136. Since the classification is statistical innature, the error in the determinations can be estimated and output as aquality indicator in each case 138. The zoning information 110 is alsoinput to this determination since it helps identify where transitionsbetween formation types might cause sudden changes in the slownesscurves which otherwise might be misinterpreted.

The track search methodology is as follows: Where Z^(k), is the set ofall measurements up through the depth k and Z(k) is the measurement atthe same depth, if at depth k n peaks are recorded, the measurement Z(k)is composed of n slownesses, n times and n semblances. Z^(k) consists ofall measurements up to and including the considered depth:

Z ^(k) ={Z(0),Z(1), . . . ,Z(k)}

If q_(i) ^(k) is the i-th global hypothesis for the measurements Z^(k),up through depth k, each q_(i) ^(k) gives one possible completeexplanation of the measurements up through the considered depth. Thusthe sequences of measurements continuing over depth (tracks) arecompletely specified and the measurements not taken into account areconsidered as false alarms. This hypothesis can be decomposed as twoparts q_(i)(k) and q_(1(i))(k). The first hypothesis, q_(i)(k),specifies at depth k, how measurement Z(k) is assigned to tracks and howtracks at depth k−1 are linked with tracks at depth k. The hypothesis,q_(1(i))(k), specifies the 1(i) hypothesis for measurement formeasurements through to depth k−1, Z^(k−1). This decomposition can berewritten as follows:

Z ^(k) =Z ^(k−1) ∪Z(k)

q ₁ ^(k) =q ₁₍₁₎ ^(k−1) ∪q _(i)(k)

This shows that it is possible recursively to generate all hypothesesfor all measurements up to the current depth, k. The tracking algorithmworks by forming hypotheses about measurements and evaluating theirprobabilities using prior probabilities model and Baye's rule. Theprobability of the global hypothesis q_(i) ^(k) based on all availablemeasurements Z^(k) can be expressed as:$P\left( {{q_{i}^{k}\left. {Z^{k}} \right)} = \frac{P\left( {{Z(k)}\left. {q_{i}^{k},Z^{k - 1}} \right){P\left( {{q_{i}(k)}\left. {{q_{l{(i)}}^{k - 1},Z^{k - 1}}} \right){P\left( \left. {q_{l{(i)}}^{k - 1}{Z^{k - 1}}} \right) \right.}} \right.}} \right.}{C_{k}}} \right.$

This equation is the general equation of the tracking algorithm withoutany simplification assumptions. The denominator C_(k) is a normalisingconstant obtained by summing the numerator of the equation over i.Conceptually, the tracking search part is completely defined by thisequation. Each time a new frame of data Z(k) is acquired, thehypothesis, q_(i)(k), at this depth is generated and appended to thehypothesis for the measurements up through k−1, the previous depth,q_(1(i)) ^(k−1) and finally the prior probability and data likelihoodmodels are used to compute P(q_(i) ^(k)|Z^(k).)

The last term of the numerator in the equation, P(q_(i) ^(k−1)|Z^(k−1))represents the probability of the parent global hypothesis which isconsidered to be known from the previous depth. The two other terms ofthe numerator equation can be obtained according to Kurien, T, 1992,Issues in the design of practical multi-target tracking algorithms.,Multitarget-multisensor tracking: applications and approaches., ArtechHouse. The term P(q_(i)(k)|q_(1(i)) ^(k−1), Z^(k−1)) is defined asP(q_(i)(k)q_(l(i))^(k − 1), Z^(k − 1)) = P(g_(r)O_(i), N_(i), Z^(k − 1))P(h_(s)O_(i), O_(l(i)))P(O_(i)O_(l(i)), Z^(k − 1))

where O_(i), the modelorder, gives the number of tracks at one depth kimplied by q_(i)(k) and O_(1(i)) is the modelorder implied by q_(1(i))^(k−1). N_(i) is the sum of all the tracks observed at the depth k andthe number of false alarm measurements implied by q_(i)(k).g_(r)(O_(i),N_(i)) is considered as the r-th way to classify N_(i)measurements as false alarms and arrivals up to O_(i) arrivals. In factthe current modelorder O_(i) and the number of measurements N_(i)determine the number of possibilities in which N_(i) measurements areclassified as false alarms and arrivals up to a maximum of O_(i)arrivals. The probability P(g_(r)|O_(i),N_(i)) is computed with theassumption that each O_(i) arrival is independent of any others and ofthe false alarms. Supposing that for one frame of O_(i) arrivals eachhas a probability p to appear at the current depth and a probability 1−pnot to appear. The false alarms are modelled as independent and theirnumber is modelled as Poisson distribution with mean 1. Considering K ofthe O_(i) arrivals present in the classification g_(r)(O_(i),N_(i)) thenthere are K arrivals and N_(i)−K false alarms and the probability is:${P\left( \left. {{g_{r}\left( {O_{i},N_{i}} \right)}{{O_{i},N_{i}}}} \right) \right.} = {{p^{K}\left( {1 - p} \right)}^{O_{i} - K}\frac{{\exp \left( {- \lambda} \right)}\lambda^{N_{i} - K}}{\left( {N_{i} - K} \right)!}}$

h_(s)(O_(i),O_(1(i))) is a function that gives the s-th possible way toassign O_(i) tracks at depth k with the O_(1(i)) tracks in the previousdepth k−1. This function gives an enumeration of the different ways tolink tracks in the current depth with those in the previous depth.Considering that all the combinations are equally probable, forψ(O_(i),O_(1(i))) possibilities, then${P\left( {{h_{s}O_{i}},O_{l{(i)}}} \right)} = {\frac{1}{\psi \left( {O_{i},O_{l{(i)}}} \right)}.}$

The classification of measurements as arrivals or false alarms,g_(r)(O_(i),N_(i)), and track assignment, h_(s),(O_(i),O_(1(i))), withthe probability models P(h_(s)|O_(i),O_(1(i))) andP(g_(r)(O_(i),N_(i))O_(i),N_(i)) give a complete determination ofP(q_(i)(k)|q_(1(i)) ^(k−1),Z^(k−1)).

The last term of the concept equation is computed as follows: Given dataZ(k) at depth k decomposed into a continuous part C(k) and a discretepart N(k): Z(k)={C(k), N(k)}.

Considering C(k) and N(k) are independent, P(N(k)|q_(i) ^(k), Z^(k−1)can be rewritten as:

P(Z(k)|q _(i) ^(k) , Z ^(k−1))=P(C(k)|q _(i) ^(k) , Z ^(k−1))P(N(k)|q_(i) ^(k) , Z ^(k−1))

where P(N(k)|q_(i) ^(k),Z^(k−1))=1 if N_(i)=N(k) and 0 otherwise.

q_(i) ^(k) completely describes the sequences of measurements whichbelong to the same tracks. To know if the sequences of measurementsenumerated as “sufficiently continuous” to qualify as real tracks ornot, the probability likelihood model for the continuous part of thedata P(C(k)|q_(i) ^(k),Z^(k−1)) must be provided.

C(k) is formed by the slowness-time co-ordinate of each peak modelled asa 2D measurement vector. Assuming that the prior probability modelconsidered is that the sequences of slowness and time are Gaussianrandom process, the correlation of slowness and time along depth ismodelled as the output of an ARMA filter with known coefficients. Tocompute this probability, Kalman filter theory is used. Scharf, L.,1991, Statistical signal processing: Wiley, is used to implement theARMA filter via a Kalman state representation and Harve), A., 1994. Timeseries models: MIT Press, to implement the Kalman update and predictionequations. Given the prior continuity model for a track I, and usingsuccessive conditioning, Kalman filter equations allow the computationof P(C^(k)|I) which is the probability of continuous measurements fromdepth 0 to k:$P\left( {{{C^{k}\left. I \right)} = {\prod\limits_{{i = 1},k}\quad {{P\left( {C(i)} \right.}C^{i - 1}}}},I} \right)$

Each of the terms in the equation can be computed using Kalmanprediction and update equations. At each depth k, Kalman predictiongives the 2D Gaussian probability distribution of the slowness-timegiven the past measurements C^(i−1). The likelihood of the data C(i)belonging to the track is computed by evaluating the Gaussiandistribution at the slowness-time point specified by C(i). Knowing thisinformation, the 2D Gaussian probability distribution of siloslowness-time at depth i+1 is predicted using C(i) and the result ofupdating of the state and covariance matrix via the Kalman updateequation. This is the final term P(Z(k)|q_(i) ^(k),Z^(k−1)) needed tocomplete the likelihood computation above.

By applying this methodology to the measurement data, measurements canbe linked together on a depth-by-depth basis to form tracks. Thosemeasurements that cannot be so linked are rejected. It is then necessaryto assign each track to a physical attribute, for example monopolecompressional arrivals, monopole shear arrivals, etc. The approach usedin this invention is to first make certain assumptions, for example thatthe first arrival track will be monpole compressional followed bymonpole shear followed by Stonely. Given these assumptions, the datamaking up the tracks is then tested to determine the probability that atrack fits these assumptions and are classified according to whichassumption shows the highest probability. The general methodology forthis determination is as follows:

D^(k) is the track data at depth k and contains the followinginformation: m, called the modelorder, that gives the number of tracksidentified at the k-th depth, the triplet (P, t, r), the slowness, timeand semblance of each track, and some linking information that provideshow a track at depth k is linked with the detected tracks at depths k−1and k+1. H_(j) ^(k) is the j-th hypothesis at depth k. The index jcorresponds to the number of hypotheses in the particular case, forexample in the case of monopole the value is 3 (for the three modes:compressional, shear and Stoneley). The posterior probability H_(j) ^(k)given the data at depth k−1 and all previous depths is defined by thefollowing equation where (.) denotes (D^(k−1), . . . D¹):$\begin{matrix}{{P\left( {{H_{j}^{k}D^{k}},\ldots \quad,D^{1}} \right)} = \quad \frac{P\left( {D^{k}\left. {H_{j}^{k},\left( . \right)} \right){P\left( {H_{j}^{k}\left. \left( . \right) \right)} \right.}} \right.}{\sum\limits_{j}\quad {P\left( {D^{k}\left. {H_{j}^{k},\left( . \right)} \right){P\left( {H_{j}^{k}\left. \left( . \right) \right)} \right.}} \right.}}} \\{= \quad \frac{P\left( {D^{k}\left. H_{j}^{k} \right){P\left( {H_{j}^{k}\left. \left( . \right) \right)} \right.}} \right.}{\sum\limits_{j}\quad {P\left( {D^{k}\left. H_{j}^{k} \right){P\left( {H_{j}^{k}\left. \left( . \right) \right)} \right.}} \right.}}}\end{matrix}$

This considers that the data likelihood for the data at the currentdepth depends on the previous data only through the hypothesis at thecurrent depth, which means P(D^(k)|H_(i) ^(k), (.))=P(D^(k)|H_(i) ^(k)).These two equations are the main equations of the classification.

The hypotheses for a given depth depend upon the number of tracks (m) atthis depth. To assign the probabilities P(H_(j) ^(k)|H_(i) ^(k−1)) somehypotheses at the current depth will be taken as consistent with thehypotheses at other depths and other will not. For example, if theclassification classifies one track at depth k−1 as a component A andclassifies a track at depth k as the same arrival A, this gives the mostprobable chance to classify the detected track as arrival A compared toany other. Therefore, P(H_(j) ^(k)|H_(i) ^(k−1))=K if the hypothesesH_(j) ^(k) and H_(i) ^(k−1) are consistent over the depth and zerootherwise. The value K is determined as follows: for M consistenthypotheses found, K=1/M. Thus K is linked with the number of consistenthypotheses found. In this case, all the hypotheses are taken intoaccount as equally probable and it is assumed that the track search parthas made no mistakes when linking measurements over the depth. Given afinite probability of error in the track search part, P(H_(j) ^(k)|H_(i)^(k−1))=p when the hypotheses are consistent over depth and 1−potherwise, with 0≦1−p≦1. In this case it is necessary to set the valueof p inside the algorithm.

The computation of the likelihood P(D^(k)|H_(j) ^(k)) assumes that thetracks are independent of each other the model used to computeP(D^(k)|H_(j) ^(k)) considers that each track is a slowness-time pointin the 2D slowness-time plane which is Gaussian distributed. The priorprobability of this slowness-time point depends upon the classificationassigned to that point by H_(j) ^(k). If the number of tracks equalsone, then the P(D^(k)|H_(j) ^(k)) is a Gaussian distribution with meanand variance given by H_(j) ^(k) and evaluated at the slowness-timepoint given by D^(k). For more than one arrival, the assumption is thatthe independency of tracks and P(D^(k)|H_(j) ^(k)) is the product of allGaussian distributions evaluated.

It will be appreciated that the tracking and classificationmethodologies described above can be applied to sonic logging dataindependently of the other processing steps described here.

The slowness logs 136 together with the appropriate error bars 138 areoutput to a finalisation process 140 which computes and selects the bestcompressional, slowness, the best shear slowness and the best Stoneleyslowness out of those available 142, e.g.:

1. Compressional Slowness:

P&S Rec., Tra., or BHC

UDP Rec., Tra., or BHC

LDP Rec., Tra., or BHC

2. Shear Slowness:

P&S Rec., Tra., or BHC

LUDP Rec., Tra., or BHC

LDP Rec., Tra., or BHC

3. Stoneley slowness:

MST Rec., Tra, or BHC

The process selects the log with the minimum error bar.

The best slowness is computed from one or more datasets using thefollowing steps:

1. A mean error is computed over the whole interval over depth for everycandidate.

2. The mean error is computed for each processing mode, i.e. Rec. orTra. If the candidate is processed by BHC mode, the mean error valuesfrom both processing modes are averaged.

3. Absent value resident in the error bar logs is regarded as a largeerror in the error computation.

4. A candidate is selected which has the minimum mean error.

5. Slowness is computed level by level based on the error bars of theselected candidate.

6. If an error at a level is bigger than a threshold, absent is assignedas a slowness at the level (no-log condition). If the candidate isprocessed by BHC mode, the procedure looks over both errors of receiverand transmitter mode. If their errors are similar, the mean slowness isassigned as the slowness, otherwise the one with the smallest error ischosen.

7. Final error bar logs associated with the selected slowness logs arederived in the same way 144.

The general multi-shot processing methodology preferred for the presentinvention is as follows:

For the case of multiple shots at depth Z_(n):

A_(i) is the waveform taken into account to process MSTC.

i is the number of a sub-array (i=1, N).

F_(i) is the ST plane computed for the sub-array i.

A frame is represented as a function of two parameters: F_(i)=F_(i)(t,S)where t is the time and S the slowness. For N sub-arrays, N=1, 5 forexample, five frames F={F₁, F₂, F₃, F₄, F₅} are obtained where F is thefamily of frames at depth Z.

The object of the multiple shot is to compute$F_{stack} = {\sum\limits_{n = 1}^{N}\quad F_{i}}$

But as there is a shift in time between two consecutive frames, it isnecessary to take into account this time shift before computing the MSTCprocessing.

Between frame F_(i) and F_(i−1), the tool moves a known distance d=Δz.For a frame F_(i)(t,S), to shift this array before computation of themultiple shot, it is necessary to apply the following operation: for thepoint A(t_(i), S_(i)) located in the ST plane, it is necessary to shiftit by:

tt _(i) =t _(i) −d·S(t _(i)).

For the totality of the frame, after time shift, this gives:

F _(i)(t,S)=F _(i)(t−(i−1)d·S(t),S).

This equation, called the time shift equation, explains how to make ashift in time between one frame and the succeeding frame. Conceptuallyto stack to ST plane of each array it is necessary to apply thiscorrection.

Starting with the equation above, it is possible to rewrite it as ageneral equation:

F _(i)(t,S)=F _(i)(t,S)*δ(t−(i−1)d·S(t),S)

where * is the convolution operator on the variable t.

To resample the frame to stack it with others, the shah function isused, defined as:${{shah}\quad (t)} \equiv {\sum\limits_{n = {- \infty}}^{+ \infty}\quad {\delta \quad {\left( {t - n} \right).}}}$

For a signal s(t) it is possible to rewrite this signal as:${s_{\tau}\quad (t)} = {s\quad (t)\quad {shah}\quad {\left( \frac{t}{\tau} \right).}}$

Thus the signal s(t) is discretised with a sampling rate τ.

Rewriting the equation of a frame shifted in time and after resamplinggives:${F_{i,\tau}^{\prime}\quad \left( {t,S} \right)} = {F_{i}^{\prime}\quad \left( {t,S} \right)\quad {shah}\quad {\left( \frac{t}{\tau} \right).}}$

Developing this equation gives: $\begin{matrix}{{F_{i,\tau}^{\prime}\quad \left( {t,S} \right)} = \quad {F_{i}^{\prime}\quad \left( {t,S} \right)\quad {shah}\quad \left( \frac{t}{\tau} \right)}} \\{= \quad {\left\lbrack {F_{i}\quad \left( {t,S} \right)*\delta \quad \left( {t - {\left( {i - 1} \right)\quad {d \cdot S}\quad (t)}} \right)} \right\rbrack \quad {shah}\quad \left( \frac{t}{\tau} \right)}} \\{= \quad {{\tau }\quad {\sum\limits_{n = {- \infty}}^{+ \infty}\quad {\left\lbrack {F_{i}\quad \left( {t,S} \right)*\delta \quad \left( {t - {\left( {i - 1} \right)\quad {d \cdot S}\quad (t)}} \right)} \right\rbrack \quad \delta \quad \left( {t - {n\quad \tau}} \right)}}}}\end{matrix}$

This defines the transformation to be applied on each frame beforestacking.

Knowing the transformation that it is necessary to apply on a frame i,it is possible to write the general equation of multi-shot processing.For a given level Z, this is:$F_{stack} = {{F_{1}\quad \left( {t,S} \right)} + {\sum\limits_{i = 2}^{N}\quad {\sum\limits_{n = {- \infty}}^{+ \infty}\quad {{{\tau }\left\lbrack {F_{i}\quad \left( {t,S} \right)*\delta \quad \left( {t - {\left( {i - 1} \right)\quad {d \cdot S}\quad (t)}} \right)} \right\rbrack}\quad \delta \quad {\left( {t - {n\quad \tau}} \right).}}}}}$

This equation is the general multi-shot equation independent of themedium (homogenous or heterogeneous). To avoid the problem of increasecomputation time due to resampling in time and time shift an alternativeis proposed that is conceptually similar to the multi-shot equation. Theidea is to include the time shift problem during the processing of MSTC.Presently, when computing MSTC, one receiver is selected as referencefor all of the processing. In the present case, the reference receiverof each sub-array is changed before computing the semblance of thissub-array. The preferred approach is to use as a reference the firstreceiver of the first sub-array, the second receiver of the secondsub-array and so on. Each time the semblance for each of thesesub-arrays is computed. This takes into account the time shift betweenthe different sub-arrays, so it is possible to stack the semblancesimmediately after computation, the time shift between the differentsub-arrays being taken into account.

After the time shift correction has been made, and the semblance foreach sub-array computed, it is necessary to stack the ST planes by usingan algebraic mean to compute the final ST plane at one considered depth.This is achieved by computing:${F\quad \left( {S,t} \right)} = {\sum\limits_{i = 1}^{N}\quad {F_{i}\quad \left( {S,t} \right)}}$

where F_(i) is the ST plane computed for one sub-array, and N is thenumber of sub-arrays used to compute the final multi-shot result.

Efficient processing of multi-shot data requires that the following bedetermined:

1. The number of sub-arrays used to compute the multi-shot.

2. Selection of good sub-arrays for computing multi-shot.

3. Automatic detection of the reference receiver for each sub-array.

For a tool with R receivers, R_(sb) is the number of receivers in asub-array. This number is selected by the user and is linked to theresolution to be obtained. By definition its lower limit is 3. Thenumber does not change the actual processing, just the physicalinterpretation of the result of the processing. N_(sb) is the maximumnumber of sub-arrays used to compute the multi-shot. This value is equalto the number of measurement depths through which the tool is moved tomake the multi-shot measurements and can be defined as:

N _(sb) =R−R _(sb)+1.

In order to know which receivers it is necessary to use to compute themulti-shot, a “geometrix matrix” is implemented. When computing amulti-shot, it is necessary to change the reference receiver for eachsub-array when computing the ST plane to take into account the timeshift correction. For N_(sb) sub-arrays of R_(sb) receivers, thegeometrix matrix is defined as: $G = \begin{pmatrix}{r_{n + R_{sb}} - 1} & {r_{n + R_{sb}} - 2} & \ldots & r_{n + 1} \\{r_{n + R_{sb}} - 2} & {r_{n + R_{sb}} - 3} & \ldots & r_{n} \\\vdots & \vdots & ⋰ & \vdots \\r_{n} & r_{n - 1} & \ldots & r_{n - N_{sb} + 1}\end{pmatrix}$

The first receiver of the first sub-array r_(n) is equal to the value ofsub-array necessary to compute the multi-shot. Therefore, r_(n)=N_(sb),the second receiver of this sub-array is r_(n+1) and so on. Each columnof the matrix can be computed as:

r _(i) =r _(n) +i, i=0, R _(sb),

Where i is the number of the receiver inside the sub-array. Thereference receiver is always the same as the first receiver of the firstsub-array. Thus if this value is computed at the beginning of thecalculation, all of the parameters for the STC processing are fixed forthis computation and it is not necessary to redefine the referencereceiver for each sub-array. Furthermore, the geometrix matrix allowseasy identification at each frame the receivers that are activated anddeactivated at that depth. Consequently, it is merely necessary tomultiply the receiver status by the current vector to know whichreceivers are to be used, thus simplifying data management. All that isneeded is to know whether or not a receiver is active for a given frame,since data from all receivers is calculated at the beginning and theresults used by circular permutation of the column of the matrix foreach new frame.

It will be appreciated that this approach can be used for any multi-shotprocessing and not just in accordance with the present invention.

The present invention finds application in the field of characterisingunderground formations surrounding borehole such as in the oil and gasindustry.

What is claimed is:
 1. A method of determining the sonic slowness of anunderground formation from sonic measurements, comprising: (i) obtainingsonic waveforms in the underground formation; (ii) determining, from thesonic waveforms, multiple values of at least one parameter related tothe sonic slowness of the formation together with an estimate of theerror in each value; and (iii) using the estimate of the error in eachvalue to select a parameter value related to slowness as representativeof the formation.
 2. A method as claimed in claim 1, wherein themultiple values of the parameter related to slowness include multiplevalues of compressional, shear and/or Stoneley slownesses.
 3. A methodas claimed in claim 1, wherein the sonic waveforms are obtained using atool having one or more transmitters and an array of receivers.
 4. Amethod as claimed in claim 3, wherein sonic waveforms are obtained intransmitter mode, receiver mode and/or borehole compensated mode.
 5. Amethod as claimed in claim 3, wherein the sonic waveforms are derivedfrom monopole or dipole measurements.
 6. A method as claimed in claim 1,wherein the step of determining the parameter related to sonic slownessincludes the step of determining a formation type parameter from thesonic waveform.
 7. A method as claimed in claim 6, wherein the formationtype parameter is derived from a compressional slowness log whereinthresholds have been applied thereby squaring the log output to fit intopredetermined bands.
 8. A method of determining the sonic slowness of anunderground formation from sonic measurements, comprising: (i) obtainingsonic waveforms in the underground formation; (ii) determining, from thesonic waveforms, multiple values of at least one parameter related tothe sonic slowness of the formation together with an estimate of theerror in each value; and (iii) using the estimate of the error in eachvalue to select a parameter value related to slowness as representativeof the formation, wherein the step of determining the parameter relatedto slowness further includes selecting processing parameters relating towell properties, formation properties and waveform processingparameters, and applying these parameters to, the processing of sonicwaveforms.
 9. A method as claimed in claim 8, wherein the processingparameters comprise at least one parameter that is selected by a userand at least one parameter that is selected or calculated on the basisof the at least one parameter selected by a user.
 10. A method asclaimed in claim 9, wherein the number of parameters selected by a useris relatively small compared to the number of parameters that areselected or calculated on the basis of the at least one parameterselected by a user.
 11. A method as claimed in claim 8, wherein thewaveforms are processed using a slowness time coherence technique.
 12. Amethod as claimed in claim 11, wherein a multishot technique is used informations in thin beds.
 13. A method of determining the sonic slownessof an underground formation from sonic measurements, comprising: (i)obtaining sonic waveforms in the underground formation; (ii)determining, from the sonic waveforms, multiple values of at least oneparameter related to the sonic slowness of the formation together withan estimate of the error in each value; and (iii) using the estimate ofthe error in each value to select a parameter value related to slownessas representative of the formation; and (iv) linking the values of theparameter related to slowness along the interval of interest into trackswhich describe the development of a particular slowness parameter inthat interval.
 14. A method as claimed in claim 13, further comprisingassigning identifiers to each track which identify a particular slownessparameter to which the track data relates.
 15. A method as claimed inclaim 14, wherein the estimate of error in each value is determined foreach track.
 16. A method as claimed in claim 15, wherein, for a givenparameter related to slowness with multiple candidates, the candidatewith the minimum mean error is selected as representative of thatparameter.
 17. A method as claimed in claim 16, further comprisingoutputting error logs for the selected candidate logs.
 18. A method asclaimed in claim 2, wherein the sonic waveforms are obtained using atool having one or more transmitters and an array of receivers.
 19. Amethod as claimed in claim 4, wherein the sonic waveforms are derivedfrom monopole or dipole measurements.
 20. A method as claimed in claim11, further comprising linking the values of the parameter related toslowness along the interval of interest into tracks which describe thedevelopment of a particular slowness parameter in that interval.